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Abstract 

Using a normal form approach described in a previous paper we derive an amplitude 
equation for a reaction-diffusion system with a Hopf bifurcation coupled to one or more 
slow real eigenmodes. The new equation is useful even for systems where the actual 
bifurcation underlying the description cannot be realized, which is typical of chemical 
systems. For a fold-Hopf bifurcation, the equation successfully handles actual chemical 
reactions where the complex Ginzburg-Landau equation fails. For a realistic chemical 
model of the Belousov-Zhabotinsky reaction, we compare solutions to the reaction- 
diffusion equation with the approximations by the complex Ginzburg-Landau equation 
and the new distributed fold-Hopf equation. 
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1 Introduction 

Near a bifurcation, the evolution of a dynamical system exhibits critical slowing down, which 
often admits a simplified description in terms of an amplitude equation. We demonstrate 
in this paper that slow modes may sometimes be described in terms of amplitude equations 
even if the underlying bifurcation cannot be realized for a given system. More specifically, 
we show that slow modes additional to the critical modes of a bifurcation may be treated in 
terms of a higher-codimension bifurcation. 

Near a Hopf bifurcation, a reaction-diffusion system may sometimes be described by a 
complex Ginzburg-Landau equation (CGLEQ), where the motion is restricted to the unfolded 
center manifold of the focus. Diffusion may tend to move the system off that slow manifold, 
but fast non-oscillatory modes may keep the state point close to it. However, if the system 
has a slow real mode, diffusion may significantly excite that mode, thus making a description 
in terms of the CGLE dubious. 

In principle, a system with a Hopf bifurcation and a slow real mode may always be treated 
by a CGLE provided we work sufficiently close to the bifurcation point. But experimental 
limitations may prevent us from operating close enough. Slow modes are likely to appear in 
complex chemical or biochemical reactions where the spectrum of characteristic times often 
spans many orders of magnitude. An example is the Belousov-Zhabotinsky (BZ) reaction 
which under certain conditions shows simple periodic oscillations. These basic oscillations 
may be described by relatively simple models like the Oregonator H|. At other operating 
points, the BZ reaction may exhibit complex oscillations (including chaos) |^, |J, which can be 
attributed to the presence of a slow real mode and cannot be described by the Oregonator. 
But the complexity of the BZ reaction can be described quite well by a four-dimensional 
extension of the Oregonator model For that more realistic model, we show in this paper 
that a description by the CGLE fails at experimentally feasible conditions. 

Treating a reaction-diffusion system by an amplitude equation saves computational re- 
sources and provides interpretations that give more insight into the character of the solution. 
It is therefore highly desirable to develop an amplitude description to replace the CGLE for 
important chemical systems like the BZ reaction or biochemical reactions such as glycolysis. 
In this paper, we shall derive an amplitude equation for reaction-diffusion systems with a 
Hopf bifurcation and a slow real mode. The equation is based on the amplitude equation 
for a fold-Hopf bifurcation even though the system need not posses such codimension-two 
bifurcation. We shall refer to the result as a distributed slow-Hopf equation (DSHE). We 
show that the four-dimensional reaction-diffusion model of the BZ reaction may be excellently 
described in terms of the appropriate DSHE. 

The development of the DSHE involves several intermediate results that are important 
in their own right. We first derive a set of amplitude equations for the class of bifurcations 
(occurring in the corresponding spatially homogeneous system) with one pair of pure imagi- 
nary eigenvalues and any number of semisimple eigenvalues zero. (Chemical and biochemical 
reactions may well have several slow modes.) 

The method used for deriving the amplitude equations builds on results previously ob- 
tained for ordinary differential equations The unfolded center manifold for the corre- 
sponding homogeneous system is parameterized by points y in the center subspace. At each 
point of (physical) space, the state point is restrained to move on the center manifold, and 
the amplitude equation is found as the equation for the modulation of the oscillations of y 
(described in a basis of critical eigenvectors) obtained from the restrained reaction-diffusion 
equation. 

1 In this paper, we use the following abbreviations: CGLE: "Complex Ginzburg-Landau equation"; DSHE: 
"Distributed slow-Hopf equation"; BZ: "Belousov-Zhabotinsky". 
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The resulting class of amplitude equations for the general case is given by Eq. (JlSj) with 
Eqs. ( p0| ) and (|2l]). From a solution to an amplitude equation, we get the physical solution 
using the transformation (E2), determined as indicated in Appendix A (see also wm). 



We treat the special case of a distributed fold-Hopf bifurcation in detail in Section 2.3 
(see Eq. ( |24| ) together with Table [l]). We note also that the CGLE is another special case; 
thus the present work contains a new derivation of that important equation. 

We modify the distributed fold Hopf equation to get an amplitude equation, the DSHE, 
applicable to systems with a Hopf bifurcation and one slow but non-critical real mode. The 
result is Eq. (|t]), which can be scaled to the dimensionless form Eq. ([n]). The precise relation 
of the DSHE to the original CGLE is required when the two equations are compared. It is 
explained in Appendix |Cj. 

The approximations in terms of various amplitude equations are evaluated for a four- 
dimensional extension of the Oregonator model mentioned above, which is described in 



Section |4.l| . The model provides a realistic description of the BZ reaction based on ex- 
tensive kinetic work |], || and describes much of the dynamical behavior of the BZ system 
qualitatively or semiquantitatively correctly (see |j| for a description of the complex dynamics 
and iQ for the bare oscillations). 

We compare the "exact" solution obtained by direct numerical integration of the reaction- 
diffusion equation with the approximations obtained as solutions to the amplitude equations, 



the CGLE on the one hand and the DSHE on the other. We find in Sections [L3| and [L4 
that whereas the simple CGLE fails, the DSHE provides a very accurate description of the 
dynamics of the reaction-diffusion system. 

The result shows the effect of a slow real mode near a Hopf bifurcation, and that a 
Ginzburg-Landau approximation cannot be considered reliable for the BZ system. The slow 
real mode is a natural (unavoidable) part of that system. A similar effect has also been 
considered by Aranson et. al Q who coupled the FitzHugh-Nagumo model || [l0| with an 
artificial slow real mode. There the slow mode interacted with an excitable system, and the 
effect was the opposite of the one reported here for the oscillatory BZ system. Random initial 
conditions normally give rise to a disordered distribution of spiral waves in the FitzHugh- 
Nagumo model, but the slow mode resulted in the development of a single spiral covering 
the entire spatial domain. 

2 Derivation of amplitude equations 

In this section, we derive amplitude equations for chemical reaction-diffusion systems for 
which the (spatially) homogeneous system is near a bifurcation at which the spectrum of 
critical eigenvalues consists of one pair of complex conjugate, imaginary eigenvalues and 
several real semisimple eigenvalues equal to zero. 

2.1 Homogeneous systems 

We consider a dynamical system, described by a vector x g K™, and depending on a set of 
parameters /x £ M s . The system is supposed to have a stationary point x s (/x) at which a 
local bifurcation occurs at fi = 0. We shall use the stationary point, x s (0), at the bifurcation 
point as origin. The homogeneous chemical system evolves according to the kinetic equation 

^ = F(x;/x)=J-x + f(x, At ). (1) 
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The Jacobian matrix at the stationary point at the bifurcation, J = DF(0, 0), has r semisim- 
ple eigenvalues and r pairs of linearly independent right and left eigenvectors u 2 and u* , 
normalized according to the biorthonormality relations 

u* • Uj = % , for i,j = 1, . . . , r. (2) 

The r-dimensional center manifold W c at = is tangent to the center subspace E° at 
x = 0, which is spanned by the r right eigenvectors Uj. Near x = 0, motion in W c (or its 
unfolding) is usually much slower than motion towards it or away from it (critical slowing 
down) , and the motion outside W c may often be considered transient and fast and therefore 
disregarded. This fact adds to the practical importance of the center manifold. 

We have previously derived simple, explicit expressions that allow us to write down the 
normal form equation for the flow on the center manifold and unfoldings of it (for fi ^ 0) as 
well as a transformation from normal form coordinates to the original physical coordinates. 
The transformation has the form of a map 

x = y + h(y,/x), where yeE c , (3) 

from £ c xR s to the (unfolded) center manifold. For oscillatory modes, the normal form allows 
a formulation in terms of "amplitudes" by a straightforward modification of the normal form 
equation. Consequently, we sometimes refer to the equation as an "amplitude equation", 
whether it is modified or not. The results of the previous work needed in the present paper 
are summarized in Appendix [A]. 

From the results of Appendix [A[ we may easily obtain an explicit expression for the 
center manifold, parameterized by coordinates in the center subspace in a basis of critical 
eigenvectors, valid for any given local bifurcation (with semisimple critical eigenvalues). We 
refer to such representation as an "amplitude transformation". It is given by Eq. (§) with 
h expanded in powers of the coordinates of y, see Eq. (44) of Appendix |A|. Observe that 



this expansion also provides an explicit expression for the vector field on W c . We shall 
make use of these two results in the next subsection to derive an amplitude equation for the 
corresponding reaction-diffusion problem. 



2.2 Extension to spatially inhomogeneous systems 

When the system described by Eq. (|l|) is modified by introducing diffusion of the involved 
physical quantities, the time evolution of the system is governed by a reaction-diffusion 
equation 

^=F(x;/x) + D.V 2 x, (4) 

where x = x(r,t) depends on a spatial position r as well as time t, and D is a diffusion 
operator, which typically, to a good approximation, can be represented by a diagonal matrix. 

In a homogeneous system, the (unfolded) center manifold is invariant under the flow. In 
chemical reaction-diffusion systems, diffusion may take the local concentration away from 
W c even if it initially was on W° everywhere. 

However, close to the bifurcation, the motion on the unfolded center manifold due to the 
chemical reactions often is much slower than motion transverse to it, as already mentioned. 
Consequently, diffusion never takes the system far away from that "slow manifold", so the 
evolution of the spatial system may be approximately described by Eq. (Q) with x restrained 
to move on the slow manifold W c (n) given by Eq. (||). 
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We may therefore simply substitute the right hand side of Eq. (||) for x in Eq. (Q) to get 
a differential equation in y, the corresponding point in the center subspace. The appropriate 
solution to the resulting equation is then transformed by Eq. (||) to the motion in the center 
manifold, the approximate solution to the reaction-diffusion equation ([I]). 

Unfortunately, direct substitution of the parameterization for the slow manifold W c (/x) in 
Eq. (f|) results in "counter-rotating terms" that prevent straightforward use of amplitudes 
for oscillatory modes. By averaging, these recalcitrant terms can be seen to be negligible close 
enough to the bifurcation point, but that method is approximate and inelegant. We shall 
therefore use a method of two times while still building on the result of the theory outlined 
in Appendix [X]. 

We now restrict the discussion to a bifurcation with a single pair of complex conjugate 
imaginary eigenvalues ± i ojq together with any number of semisimple zero eigenvalues. There 
are r such critical eigenvalues (counting multiplicity) and r linearly independent right eigen- 
vectors and corresponding left eigenvectors satisfying the biorthonormality relations (||). The 
period of the oscillations near x = 0at/x = 0isT = 2n /uj . 

We express the solution y(r, t) as the modulation of harmonic oscillations of frequency ujo 
for the oscillatory mode and seek an equation for the modulation. To that end, it is convenient 
to view y as a function of two independent time variables, r and 9, with 9 accounting for the 
harmonic oscillation (of the oscillatory degrees of freedom) and r describing the modulation. 

We define 

r 

y(r,T,0) = e je -z(r,r) =^e^(r,r)u 4 , (5) 

i=l 

where r and 9 are given functions of time. At the end of the calculation, we shall choose 

r = 9 = t. (6) 
To obtain the correct equation for y(r, r, 6), we use Eq. (g) initially as 

dy = dyd£ dy dr = dy dy 

dt 06 dt dr dt 86 dr' { ' 

but otherwise we use the variables r, r, and 9 consistently throughout (so no confusion should 
arise from using the symbol y for the two different functional forms). 

We now view x as a function of r, t, and 9 through Eq. (^), and using Eq. (|?|), the 
reaction-diffusion equation may then be written as 

^ = (J.x-g) + f(x, / x) + D.V 2 x, (8) 

with 

x(r, r, 9) = y(r, r, 9) + h(y(r, r, 9), /x). (9) 

For a given physical (or chemical) problem, the transformation h is a well defined function 
of y and /x. So for a given \i, the second term of Eq. (||) and its derivatives with respect to 
r, t, and 9 are completely determined by the function y(r, r, 9). Consequently, these terms 
do not contain any "new" information about the dynamics not already contained in y(r, r, 9) 
and the transformation h(y, fj,). 

These remarks suggest that we try to eliminate as much as possible of the nonlinear term 

of 
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in Eq. (|8[) . We observe that 



i—l 

implying that we may extract an equation for by multiplying Eq. (Js[) by e~ Xi9 u* (from 
the left) and averaging the result over a period T of 8. The left-hand side of Eq. (||) then 
simply gives ^f, since 

We prove this property in Appendix ||. There we also show that 



— / c 



«J-x-— )^ = 0, (13) 

implying that the two terms in parentheses disappear from Eq. (||) by the above operation. 

As argued in Appendix |A|, we may use a power series expansions of x and f (x, /z) for x 
on the slow manifold W c (fi). Here we first note that 

yP = JJyf = CX P £>A^)IK ! = ° PA9zP (14) 

in concise notation. 

For f (x, /a), we get using the expansion (|is| ) 

i / T c- A ^u*-f(y + h(y^))^-^E u *- f p q z > q [ T c {p X - x > )e de 

1 J 1 pq JO 

= x;v^M q , (is) 

pq 

in which / pq = u*-f pq is the i'th component of the vector coefficient of the expansion of 
f(x, /i) on W c (n) whereas the sum is taken over all sets (p,q) for which the resonance 
condition for the i'th component 

P' A = E p i A i =Ai (I 6 ) 



is satisfied. All other terms vanish because of the integral of the exponential in Eq. 
Note that critical eigenvalues are either pure imaginary or zero. For clarity, the summation 
over the resonant terms for the i'th component is marked with a superscript i in the sum in 
Eq. ( |l5| ) . For the diffusion term of Eq. (g) , we similarly get (using the expansion (Q) from 
Appendix |Aj) 



^ V A * V*'D-V 2 (y + h(y, M ))^ 



pq " pq 



E A,'V 2 ^. + E l <-D-h P qM q V 2 zP. (17) 



pq 
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Here 

!>,, ",Du. (18) 

and the sum over j is taken over all terms for which Aj = Ai, and the sum over p and 
q includes all resonant terms for the ?'th component (and no other). By using Eq. (|T^), 
Eq. (|l3|), Eq. (|l5|), and Eq. (|l7|), we obtain an amplitude equation for the modulation as a 
set of coupled equations for the set of all coefficients Zj(r, r) defined in Eq. (0). 



at 



^f pq zP^ + ^^V 2 z J+ ^P«V 2 zP. (19) 



p 



Here 




and 



E|q|>0 h 5 J qA tq ) J A J f = A » 

Aj 7^ A.; 



p " \o, | P | < i, (21) 

in which |p| = J2jPj- As before, a superscript i restricts a summation to resonant terms 
for the z'th component. In Eq. (po|), the index set Sj has components described by the usual 
Kronecker delta, (Sj)k = Sjk- 

In Eq. ([l9|), we have identified r with the "real" time, t, since we no longer need to 
distinguish between the two "formal" time scales. From a solution Zi(r, t), i = 1, . . . , r, to 
Eq. ( |l9|) we get an (approximate) solution to the reaction-diffusion equation (^) as 

x(r,i)=y(r,i)+h(y(r,i),/*), (22a) 

r 

y(r,t) =^z,(r,t)e A - t u l (22b) 

i=i 

Here we note that the transformation h produces anharmonic terms in 9 in the plane of 
oscillations as well as components off the center subspace E c . Of course, these terms are also 
modulated (by products of powers of the amplitudes Zi) and are therefore not periodic in t. 

For a spatially homogeneous system, the diffusion terms disappear from the amplitude 
equation (|lS|), and the result therefore agrees with the general expression for a homogeneous 
system as quoted in Appendix [A| In general, the amplitude equation contains diffusion terms 
that modify the motion on the slow manifold. The sum over j in Eq. ( |l9|) represents the linear 
diffusion terms which include couplings between bifurcating modes corresponding to the same 
eigenvalue. But note in particular that the oscillatory modes (with eigenvalues ±iwo) do not 
couple to any of the real modes — only the real modes may couple to each other. 

The nonlinear diffusion terms, the last sum in Eq. (|l9|), arise because the motion is re- 
strained to the slow manifold; diffusion is linear in x but we describe the dynamics in terms 
of points y£E c , and x is clearly not linear in y. It is to be expected, that the linear diffusion 
terms will dominate over the nonlinear diffusion, and these higher order terms may therefore 
be neglected in first approximation. Below we show examples which are excellently described 
using linear diffusion alone. In some cases, the nonlinear diffusion may be important, how- 
ever. For example, this may be the case in situations where the real part of any of the 



(20) 
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linear diffusion coefficients becomes negative — a situation which easily arises in biochemical 
systems where the diffusion constants may differ by several orders of magnitude. 

In Section |2.3| , we shall exhibit the amplitude equation for a fold-Hopf bifurcation (which 
has one oscillatory and one real mode) as a special case of Eq. (JTs|) . We also note how the 
well-known complex Ginzburg-Landau equation (CGLE) appears as a special case of the 



result ( |19[) . Thus the present work contains a novel derivation of the CGLE. 

The fold-Hopf bifurcation will be used to solve the reaction-diffusion equation for a re- 
alistic model of the Belousov-Zhabotinsky reaction — we shall compare the solution obtained 
from the new equation with those of the reaction-diffusion equation itself and of the corre- 
sponding CGLE. In addition, this example also serves as an illustration of the practical use 
of Eq. ©. 

2.3 Distributed Fold-Hopf equation 

In the following, we shall focus on the special case of Eq. (pj|), where there is just a single 
eigenvalue in addition to the imaginary pair. This bifurcation is called a fold-Hopf bifurcation 
corresponding to the case where a Hopf bifurcation coincides with a single non-degenerate 
real bifurcation. The linear diffusion terms contain no coupling between modes, so the sum 
over j in Eq. ( |l9|) reduces to just one term, du\7 2 Zi with 

d u =u*-D-(u J +^h 5iqM t i). (23) 
q 

Here we shall exhibit the amplitude equation ( |l9| ) to lowest non-trivial order in fi (i.e. in- 
cluding only terms essential to the unfolding of local terms) implying that the diffusion term 
du reduces to du = u*D-Ui. 

To lowest order, we need only terms linear in fi in general, so we get similar independent 
terms from each component of /x. Therefore, we shall exhibit the result for a single scalar fi 
even though one is often interested in using two independent parameters for bifurcations of 
codimension two: the generalization is straightforward. 

For convenience, we simplify the notation by assigning Ai,2 = ±i^o and A3 = for the 
eigenvalues at the bifurcation point, and Ui = u, U2 = u, and 113 = v for the corresponding 
right eigenvectors. In this case, we need only consider the amplitudes for w = Z\ and z = Z3 
since complex conjugation gives Z2 = w. 

If we include all terms up to third order in w and z for q = and only terms essential 
to the unfolding for q > and neglect nonlinear diffusion, we then obtain the following 
amplitude equation for a fold-Hopf bifurcation in a reaction-diffusion equation ([|) 

w = (Ti/iw + gowz + gi\w\ 2 w + g^wz 2 + d w \7 2 w, (24a) 
i = po/jL + c \w\ 2 + ciz 2 + c 2 \w\ 2 z + c 3 z 3 + d z V 2 z. (24b) 

Here we have introduced the following compact notation 





= u* 


D u, 


d. 


= V* 


D v 




= u* 


flOOli 


Po 


= V* 


foooi 


9o 


= u* 


fioio, 


9i 


= u* 


f2100 


CO 


= V* 


fllOOj 


Cl 


= V* 


f0020 



g 2 = u*-fio2o, (25) 

C2 = V -flllO, C3 = V -fo030- 

The results derived for the fold-Hopf bifurcation are summarized in Table [l| 
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The amplitude equation for a reaction-diffusion system that undergoes a simple Hopf 
bifurcation is easily obtained from the general result Eq. ( |l9| ) following the same procedure. 
As a result, we obtain the well-known complex Ginzburg-Landau equation (CGLE) 

w = ai/iw + g\iu\ 2 w + dV 2 w. (26) 

The details of the results for the CGLE and explicit expressions for calculations of the 



coefficients in Eq. (|26J) are summarized in Table |2j. The parameter d in Eq. (26) is identical 
with d w , shown in Eq. (p3). 



3 Handling a slow noncritical mode 

An amplitude equation very similar to the distributed fold-Hopf equation ( f24|) can be derived 
for a system with a Hopf bifurcation and a slow real mode, even if only a Hopf bifurcation 
can be realized (near the desired operating point or anywhere at all). The result may be 
viewed as an extension of the complex Ginzburg-Landau equation. 

3.1 The distributed slow-Hopf equation 

We consider a system depending on a scalar parameter jx, having a Hopf bifurcation (for the 
homogeneous system) at /x = 0. The critical eigenvalues are ±iwo and the corresponding 
eigenvectors are u and u. All other eigenvalues have negative real parts. One real eigenvalue 
is small, |Ao| <C loq (eigenvector v) whereas all others satisfy — Rc{A} luq. 

We may describe the system with a partial unfolding of the fold-Hopf equation, translated 
so that the origin is at the stationary point at fi = 0. Unfolding of the Hopf bifurcation is 
governed explicitly by /i and needs no further discussion. The unfolding and translation from 
the fold-Hopf bifurcation to the particular Hopf bifurcation studied, can be accounted for by 
evaluating the fold-Hopf parameters at the Hopf bifurcation point, (x, = (0, 0) instead of 
at the (fictitious) fold-Hopf point. In addition, the term poA 4 °f Eq. ( |l9| ) disappears because 
of the translation to the new stationary point, and a term \qz appears explicitly, which was 
absent at the fold-Hopf bifurcation because Ao vanishes there. The relation of unfoldings to 
parameters of an amplitude equation and amplitude transformation, evaluated at a stationary 
point of an unfolded equation, has been discussed in section VII of [[|. 

The result of the modification is an equation of exactly the same form as Eq. ( |l9| ) except 
that \$z replaces pofi. The coefficients are given by exactly the same formulas, Table EL 
except that the Jacobian matrix and the other derivatives of the vector field are evaluated at 
the Hopf bifurcation point, and so are the eigenvalues and eigenvectors of J. We shall refer 
to the equation as a distributed slow-Hopf equation (DSHE). Keeping the most important 
(lower order) terms only, the DSHE takes the form 

w = ai/iw + gowz + gi\w\ 2 w + d w \? 2 w, (27a) 
i = Aoz + coM 2 + d z V 2 z, (27b) 

which we use in Section [J to describe a realistic model of the BZ reaction-diffusion system. 

The DSHE @ is a generalization of the CGLE (||). It contains the simple CGLE as a 
special case when |Ao| 3> Re{<7i}/i, as we show in Appendix |^: there we discuss the relation 
between the two equations, and in particular we derive an explicit relation between the 
parameters g\ and g. These parameters are in general quite different despite their deceptively 
similar definitions. The two amplitude equations are evaluated and compared numerically 
with each other and with the underlying reaction-diffusion equation in Section H. 



Reaction-Diffusion Systems with Slow Real Modes 



10 



3.2 Scalings 

Near a bifurcation, an evolution equation often scales with the bifurcation parameter(s) to 
lowest order (s) in a simple way. This property admits a reduction of the (truncated ) equation 
to a form independent of the distance(s) from the bifurcation point. Examples of reduced 
equations are the Ginzburg-Landau equation Jl^|, the Kuramoto-Sivashinsky equation ||l3| , 
(li| , and the Swift-Hohenberg equation Jjl|. 

For the fold-Hopf bifurcation, a complete scaling requires two carefully chosen bifurcation 
parameters. The DSHE ( gjj ) cannot be scaled to a form independent of the bifurcation 
parameter. Nevertheless, it is useful to introduce a partial rescaling of Eq. ( p7| ) with the 
substitutions 

w = yfjiw' , z = (j,z', t = t'/fj,, r = r'/ v //Z (28) 

to obtain an equation in the new, primed variables. If we drop the primes for simplicity, the 
scaled equation becomes 

it; = <j\w + gowz + gi\w\ 2 w + d w \7 2 w, (29a) 
Hz = X z + co\w\ 2 + u.d z \7 2 z. (29b) 



For analysis of solutions and numerical computations it is convenient to reduce Eq. ( |29j ) 
to dimensionless form by further scalings through the substitutions 



-u/exp(i^'), z=°-\z\ t = t'/a[, r = J %r', (30) 
V Si \ So V °"i 

in which superscripts r and i denote real and imaginary parts respectively. The result is an 
equation in the primed variables. If we skip the primes for convenience, the resulting rescaled 
equation takes the form 

w = w + (1 + i-f)wz - (1 + ia s )w\w\ 2 + (1 + i/3)V 2 w, (31a) 
ei = \ z + k\w\ 2 + t5V 2 z, (31b) 

where 

Si a w So Si a w 

The scalings used here are similar to those traditionally used for the CGLE reducing it 
to the form 

w = w - (1 + ia)w\w\ 2 + (1 + i/3)V 2 w. (33) 



Here j3 is given by Eq. ( p2| ) and a = The parameters a and a s differ because g and g\ 
do, as discussed further in Appendix [C| 



4 Comparison of solutions 

The first successful attempt to model the oscillatory behavior of the BZ reaction was the 
Oregonator model suggested by Field and Noyes M. It exhibits both sinusoidal and relax- 
ational oscillations in semi-quantitative agreement with many experiments. However, the 
Oregonator model does not reproduce more complex dynamics such as quasiperiodic oscilla- 
tions or chaos. 
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4.1 Model of the BZ reaction 

To incorporate complex phenomena of the BZ reaction into the modeling many different vari- 
ants of the Oregonator have been considered. Here we use the modified Oregonator recently 
suggested by Wang et. al || to explain very complicated transient phenomena observed in 
the BZ reaction. The model, which we refer to as the 4D Oregonator, is based upon the 
following chemical scheme 

BrOg + Br" + 2H+ 

HBr0 2 + Br" + H+ 

BrO^ + HBr0 2 + H+ 

2HBr0 2 

BrMA + Cc 4 + 

MA + Cc 4 + 

BrMA 

Note that the first six reactions are equivalent with the reactions in the Oregonator model, ex- 
cept for the fact that bromomalonic acid (BrMA) appears explicitly in the scheme. The mech- 
anistic assumption is that HBrO reacts immediately with malonic acid to produce BrMA. In 
order to explain the complexity observed in the closed system it is essential to add a reaction 
which corresponds to a removal of BrMA without a simultaneous production of Br - . This 
additional feature is included in reaction 7 of Eq. 

Now, if we regard the concentration of BrMA as an additional dynamical variable and 
assume non-homogeneous spatial conditions, the above set of reactions gives rise to the 
following four-dimensional differential equation 

k x AB 2 Y - k 2 HXY + k 3 AHX - 2k 4 X 2 + D x VjX, 
-k ± AH 2 Y - k 2 HXY + k 5 UZ + D Y y 2 Y, 

(35) 

2k 3 AHX - k 5 UZ - k 6 MZ + D z V 2 r Z, 

k x AH 2 Y + 2k 2 HXY + k 4 X 2 - k?_,UZ - k 7 U + DuV 2 U, 

where the following short-hand notation is introduced for the concentration of the chemical 
species: X = [HBr0 2 ], Y = [Br"], Z = [Cc 4 +], U = [BrMA], A = [BrOg], H = [H+], and 
M = [MA]. As in the Oregonator model, we regard A as a constant in order to study the 
system under open conditions. Besides having a supercritical Hopf bifurcation with small 
sinusoidal oscillations, the model also contains an additional slow time scale which primarily 
is related to the kinetics of BrMA. The values of the various constants in Eq. ( |35|) are shown 
in Table || 

4.2 Bifurcations and stability of waves 

For the comparisons, we need to select a suitable operating point for the 4D Oregonator, close 
to a supercritical Hopf bifurcation. As parameters, we have chosen the fixed concentrations 
[BrOJT] and [H + ]. All other parameters have been fixed using the values shown in Table ||[ 
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The locus of Hopf bifurcations in the plane of the free parameters are shown in Fig. [l]a. 
The bifurcation diagram consists of two apparently separate curves, which are connected for 
a very large and physically unrealistic value of [BrO^], however. Both branches describe 
supercritical Hopf bifurcations. The filled circle on the upper branch indicates a Hopf bi- 
furcation point which has been chosen as reference point viz. [BrOTjnopf = 0.3662 M and 
[H+Jnopf = 1.3416M. This reference point is characterized in Table |j. Since it is necessary 
to work at a non-zero distance from the bifurcation point, the actual operating point will be 
chosen close to this reference point. The variation of the Ginzburg-Landau parameters a and 
(3, and the functional dependence (3(a) are shown in Figs. 0b-0d. In all three cases, numbers 
©-® identify curves with the branches in the bifurcation diagram. 

The reference point shown with a filled circle in Fig. [l] has been chosen with a view to 
the stability properties of waves. Reaction-diffusion equations may support several types of 
waves. Near a supercritical Hopf bifurcation, in particular, there may exist stable spiral waves 
as well as spatio-temporal chaos (chemical turbulence), depending on parameter conditions. 
Usually the CGLE approximation faithfully reproduce such behavior |p^| , and it is convenient 
to relate wave stability to the parameters, a and /3, of the CGLE. 

For the CGLE, the stability of spiral waves and their transitions to turbulent solutions in 
the a, /3-parameter plane have been studied numerically and theoretically in Jl7|, |l8, 19 1. One 



key observation in these studies is that the transition to turbulence in the CGLE is closely 
related to the concepts of convective and absolute instability. For spiral wave solutions, it 
was found that the boundary of the region for absolute instability is very close to the curve 
where spiral wave solutions exhibits a transition to turbulence. In addition, the region which 
separates the onset of convective instability of the spiral waves, to some extent, may be 
approximated by the Eckhaus instability border [ pTs| , [Ts| |. The parameter space may therefore 
be divided into several subregions described by convective and absolute instabilities as in 
Fig. [j]d. In the next section, we shall use these concepts in a discussion of suitable operating 
points to be selected for a comparison of the various equations, set up for the 4D Oregonator. 

4.3 CGLE versus the reaction-diffusion equation 

We choose the operating point well outside the region of absolute instability (AI in Fig. [j]a) 
where we expect to observe stable spiral waves for the CGLE. Observe, that the chosen point 
lies close to the region of convective instability (EH), but this has no importance to the 
local stability properties of spiral waves. The Ginzburg-Landau parameters evaluated at the 
reference point are shown in Table 0. 

The 4D Oregonator must be run at a non-zero distance from the Hopf bifurcation point. 
We choose the parameter values [BrO^"] = 0.3625 M and [H + ] = [H+Juopf- This corresponds 
to an amplitude of the uniform oscillations, which is 7.5% of the stationary Ce 4+ concentra- 
tion. We emphasize that the CGLE is scaled and hence independent of the distance from the 
bifurcation point but the reverse scaling together with the amplitude transformation account 
for the distance actually used. 

The most systematic way of comparing the solutions to the amplitude equation and the 
reaction-diffusion equation is to transform the space and time depend ent c omplex amplitudes 
back to concentration space using the scalings described in Section |3^ and the amplitude 
transformation determined in Table 0. The transformation is nonlinear, and the nonlinear 
terms are sometimes important po|| . 

For the present study, the linear part of the transformation 

x(r,t;/i) =x s (/i) + u;(r,t;^)e i " ot u + c.c. , (36) 
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is sufficiently accurate. Equation ( p6[ ) is expressed in terms of unsealed variables w, r, and 
t. In addition, the unfolding terms describing the variation of the stationary point x s (/i) 
with the parameter fi are also incorporated. As an alternative to transforming amplitudes 
to physical variables for the comparison, one may calculate the amplitudes from the physical 
variables (concentrations) obtained from the reaction-diffusion equation. This is easy in the 
linear approximation to the amplitude transformation: 

w(r,t;n) = u*-[x(r,t;/i)-x,0')]e- iw °'. (37) 

The representation (p^) is very convenient since the magnitude \w\ of the (unsealed) amplitude 
forms a very clear representation of an oscillatory wave. 

In the comparisons, Figs. || and [| we exhibit the magnitude \w(r, t; fi)\ as well as selected 
components of x(r, t; fi) (the concentrations of Ce 4+ and BrMA) for solutions to the reaction- 
diffusion equation (RDE) and the CGLE respectively, using initial conditions suggested by 
Kuramoto (we choose \w(r,t; fx)\ proportional to |r| and arg(u;(r, t; /J,)) = arg(r), the 
polar angle). 

As is clear from Fig. ||, the CGLE exhibits the behavior expected from the diagram 
Fig. [j]d, a spiral wave solution evolving from the center of the grid. However, this very 
regular solution shows virtually no resemblance with that of the reaction-diffusion equation, 
Fig. ||. Focusing on [Ce 4+ ], we see that a small spiral is initially formed at the center of 
the grid. After a few windings, the spiral is no longer able to maintain its structure, and a 
turbulent pattern develops. This is especially clear in the amplitude plot where a "shock wave" 
is emitted from the spiral core. When the shock hits the boundary of the domain, many small 
spirals are generated, the symmetry of the initial state breaks, and fully developed turbulence 
is reached. Here numerous small spirals are constantly created and annihilated, resulting in 
a disordered motion typical of chemical turbulence, which also has been observed in certain 
parameter regions for a reaction-diffusion equation based on the simple three dimensional 
Oregonator model |16|. 



4.4 DSHE versus the reaction-diffusion equation 

The failure of the CGLE is clearly a result of the presence of a slow real mode in the 
4D Oregonator, which violates a condition underlying a derivation of the equation: all motion 
towards the plane of the oscillations must be fast compared with the motion in the plane. 
We shall therefore try to describe the reaction-diffusion system with a distributed slow-Hopf 
equation, DSHE, which takes the slow real mode explicitly into account. 

We work at the same operating point as used in the previous section so the results can 
be immediately compared with those for the CGLE and the reaction-diffusion equation. The 
parameters of the DSHE are evaluated at the Hopf bifurcation point using the formulas of 
Table |. They are shown in Table Note that the parameter g\ is quite different from the 
corresponding CGLE coefficient g. This important difference is discussed in Appendix ^| 

We want to compare the solution to the DSHE with that of the reaction-diffusion equation 
first of all. For the latter, we might use Fig. ||, but to illuminate the role of the slow real 
mode, we shall calculate the real amplitude z as well as the concentrations of Ce 4+ and 
BrMA. 

The concentrations are obtained from the scaled amplitudes by reverse scaling and trans- 
formation back to concentration space with the amplitude transformation described in Table [j]. 
As before, we consider only the linear part of the transformation here: 



x(r, t; /i) = Xaifi) + (w(r, t; ^)e iuJot u + c.c.) + z(r, t; fj,)v, 



(38) 
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(shown in terms of unsealed variables). In the linear approximation (|38|) , we may alternatively 
calculate the amplitudes from a solution to the reaction-diffusion equation as 



w{r, t) = u* • x(r, t; fx) - x s (/i) 



- i ujq± 



z{v,t) 



x(r,t;/i) - x s (» 



(39a) 
(39b) 



The linear approximation (39a) to the complex amplitude w is identical to the one calculated 
for the CGLE, Eq. ( |37| ) (but the two amplitudes differ if nonlinear terms of the amplitude 
transformation are taken into account). Its magnitude \w\ has already been shown in Fig. ||, 
but to save space, we shall not exhibit \w\ for the DSHE. Instead, we show the real amplitude 
z together with [Cc 4+ ] and [BrMA] in Fig. |] for the reaction-diffusion equation and Fig. || 
for the DSHE. Initial conditions for the simulations are the same as in the previous section 

@- 

The agreement between the two solutions is striking. Even quite subtle details agree 
for a long time although, eventually, the details of the turbulent patterns must disagree. 
However, the general character of the solutions continue to agree. Thus, the improvement of 
the DSHE-approximation over that of the CGLE, Fig. [| is evident. 

In conclusion, it is possible to treat a reaction-diffusion problem near a Hopf bifurcation 
using an efficient amplitude description, even for systems like the BZ-reaction that contain 
a very slow real mode as an intrinsic part of the chemistry. Essentially, the DSHE treats the 
Hopf bifurcation (of codimension one) as an unfolded fold-Hopf bifurcation of codimension 
two. 



4.5 Role of the slow real mode 

A conspicuous feature of the solutions to the reaction-diffusion equation, Fig. [|, or the DSHE, 
Fig. |], is the similarity of the patterns for [BrMA] and the amplitude z. In fact, the two 
patterns are almost identical either in Fig. [| or in Fig. [|. 

This similarity can be explained in terms of the eigenvectors of the Jacobian matrix 
whose numerical values, evaluated at the stationary point at the Hopf bifurcation, are shown 
in Table |^. We note that the eigenvector v for the slow real mode is directed almost exactly 
along the [BrMA] -axis. 

Another noticeable feature of the solutions in Figs. || and || is the different characters of 
the components [BrMA] and [Ce 4+ [. Clearly, the spatio-temporal evolution of the BrMA pat- 
tern resembles that of \w\ much more (see Fig. |J). Again, the eigenvectors, Table [j], can throw 
some light on this feature. In the linear approximation, the left eigenvectors u* and u* (or 
the real and imaginary parts of u*) define the projection of any vector onto the plane of 
oscillation. From Table [j] we see that the [BrMA] "component" of u* is very small compared 
with the other components of u* . This implies that [BrMA] virtually does not participate in 
the basic oscillations of frequency uq. Thus, apart from scale, the amplitude z can almost be 
identified with [BrMA] - [BrMA] . 

To understand the behavior of [BrMA] and the role of the amplitude z, we first watch the 
local variations of Re w and z at an interior grid point from the simulation shown in Fig. ||. 
We observe that the amplitude R of the oscillations Re w appear to vary in roughly the same 
way as z. 

Setting w = i?exp(i 0), gives the following differential equation for R and z in the homo- 
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geneous DSHE 

R= (a[+ g r a z)R- g[R 3 , (40a) 



fiz = X Q z + c R 2 , (40b) 



where superscript r denotes the real part. To put Eq. ( |40[ ) into dimensionless form we intro- 
duce the variable change 

R=J-^R', z=^z', and t=\t'. (41) 
V a\ 9l v\ 

Skipping the primes for simplicity then gives the following dimensionless equation 

R = R + Rz-R 3 , (42a) 
. X z + kR 2 g T 

z = , where k = — Co — . (42b) 

9i 

Apart from (0,0), this scaled differential equation admits the non-trivial stationary solu- 
tion 

;? s = (i + fP, * = -7xr. (43 ) 

which corresponds to a uniform oscillatory solution to the reaction-diffusion equation. A 
numerical solution to Eq. (Q) is shown in Fig. (t| Here, we observe that the point (R s ,z s ) 
corresponds to a stable focus in the (R, z)-plane — the trajectory exhibits a spiraling motion 
converging to the stationary state in agreement with the oscillatory modulation of the am- 
plitude observed in Fig. ||. If we include diffusion into the description, this will serve as a 
constant perturbation, which keeps the solution away from the stationary solution (i? s ,z s ) 
and thereby preserves the modulation of the amplitude. The interaction between the ampli- 
tudes R and z therefore introduces an additional frequency into the system, as observed in 
the torus- like behavior of the time series for Re w, and destroys the stability of spirals. 



5 Conclusion 

We have derived an amplitude equation that provides an efficient and accurate, approximate 
description of a reaction-diffusion system near a fold-Hopf bifurcation, as well as general- 
izations to systems with several critical real modes. The derivation is based on a normal 
form transformation combined with the use of two distinct time variables to handle slow 
modulations of oscillations. It does not rely on different time scales. 

The derivation exhibits a simple geometrical interpretation of the amplitude equation 
in terms of the center manifold and center subspace, and it provides an explicit relation, 
the amplitude transformation, between amplitudes and physical variables. An amplitude 
equation is obtained as a truncation of a form containing infinite series, and there is a choice 
of including more or less terms. Similarly one may include more or less terms in the amplitude 
transformation. Observe that all coefficients can be calculated from the results of the theory. 

The main result of the paper is the general amplitude equation Eq. ( |l9[ ) for a Hopf bifur- 
cation and any number of "semisimple" critical real modes. It contains all the higher order 
terms mentioned, and in practice, it must be suitably truncated for any particular use. We 
note in particular that Eq. Ji~g| ) includes nonlinear diffusion terms. 
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The special case of a fold-Hopf bifurcation is treated more explicitly in Section 2.3. Here 



the amplitude equation (24) is a selected truncation of Eq. ( J19[ ) for the case of a single real 
critical mode. The coefficients ( p5| ) are given explicitly by Table [j] in terms of coefficients h pq 
obtained by solving the linear equations shown. 

The fold-Hopf bifurcation has codimension two, and we would not expect to find a real- 
ization of this bifurcation in any particular reaction-diffusion system. Nevertheless, we show 
that a description based on the dynamics of a fold-Hopf bifurcation is much more versatile 
than the codimension initially suggests. 

We show that the machinery used to derive the distributed fold-Hopf equation ( |24l ) can 
be used to derive what we have called the distributed slow-Hopf equation d^), the DSHE, 
applicable to a reaction-diffusion system with a simple Hopf bifurcation together with a slow 
(near-critical) real mode. The DSHE (|27]) differs from the general amplitude equation (|l9| ) 
by including the linear term Xqz which is absent at criticality. In addition, all its coefficients 
are evaluated at the Hopf bifurcation. It includes terms sufficient for the particular reaction- 
diffusion problem discussed: a four-dimensional model of the Belousov-Zhabotinsky (BZ) 
reaction exhibited in Section 11 , the 4D Oregonator. 

The purpose of looking at the 4D Oregonator ( |35| ) is first of all to test the DSHE by 
comparing solutions to it with solutions to the reaction-diffusion equation and to the corre- 
sponding CGLE, all computed at the same operating point of Eq. d3S|). 

The conclusion from the comparison (Figs. [| and ||) is that the DSHE describes the 
reaction-diffusion problem almost quantitatively in the present case. In contrast, the CGLE 
fails completely (Figs. [2] and ^). So using the DSHE, one may reduce the numerical work to 
a small fraction (typically of the order 1:100) of that required for the corresponding direct 
integration of the reaction-diffusion equation. For the operating point considered, the behav- 
ior of the CGLE agrees with the fact that the coefficients a and f3 of the equation lie in a 
region where plane wave solutions are locally stable to long wave perturbations. This spuri- 
ous stability is a consequence of the neglect of all dynamics out of the plane of oscillations. 
The correct behavior, as represented by the DSHE, can be explained by a stability analysis 
of plane wave solutions to the DSHE ( p7| ) which we have carried out using a modification of 
the approach described in Jl^ |. We shall report these results in a separate paper. 

Another reason for treating the particular reaction-diffusion system ( p5| ) , in terms of the 
DSHE, Eq. (|27|), is the importance of the BZ-reaction, which it models. The BZ-reaction 
is perhaps the most studied nonlinear chemical reaction, and happens to have a very slow 
near-critical real mode almost everywhere in parameter space. (In fact, we have not found 
a fold-Hopf bifurcation anywhere in the model.) Furthermore, the slow real mode is very 
important to the dynamics, being associated with complex oscillations and chaos. So for a 
realistic treatment of a BZ-based reaction-diffusion system, the slow real mode cannot be 
neglected. 

Recently, qualitative comparisons between solutions to the CGLE and experimentally ob- 
served spiral waves in the ferroin catalyzed BZ reaction were reported by Ouyang and Fles- 
selles [ElJ. The ferroin catalyzed BZ reaction is not so well understood and characterized as 
the Ce~catalyzed version. Although a large part of the reactions involving bromomalonic 
acid are common for the two versions of the BZ system, reactions like the fifth of Eq. J34] ) 
do depend on the catalyst. Nevertheless, the possibility of slow real modes should be kept in 
mind, and the applicability of the CGLE cannot be taken for granted. 

An advantage of the CGLE is that its coefficients a and j3 can be obtained directly 
from experiments [^2| although the methods still need some refinement for higher accuracy. 
Unfortunately, measuring the coefficients of Eq. ( p7| ) is much more difficult. We may possibly 
get some progress using more general perturbation methods based on control theory, which 
we presently are developing. Otherwise one must work through optimized models of the 
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reaction considered. 

In a case where several real modes are close to criticality, one must derive the appropriate 
amplitude equation from the general expression (|l9|) . This means calculating explicit expres- 
sions for the necessary coefficients as in Table [l]. At the same time the linear equations for 
the coefficients of the amplitude transformation must be solved to the order desired. These 
steps are quite straightforward using software packages capable of performing symbolic ma- 
nipulations such as Mathematica or Maple |^3j |24|| . Even though the multitude of coefficients 
to be calculated and equations to be solved grows rapidly with the number of slow modes, 
the gain in speed by using an amplitude equations may well make a solution feasible that in 
practice is impossible using the reaction-diffusion equation directly. 



A Amplitude transformation 

In [U, we derived an expression for the unfolded center manifold W c (fi) in terms of coordi- 
nates yi of y <G E c in a basis of critical eigenvectors Uj , namely 

r 

x = y + h(y,/x) = ^y 4 u l +^h pq y p A t q (44) 

i=l pq 

in which y p = Yl t vT an d M q = EL l^k ■ The coefficient vectors h pq are solutions to the 
linear equations 

r 

( J - ' h pq = Q P * Pq , (45a) 

R p h pq = 0. (45b) 

in which I is the unit tensor, R p projects onto the resonant subspace and Q p removes the 
resonant components of any vector x: 

R p .x=(I-Q p ).x = 2«.x)u i . (46) 

i 

The sum in Eq. ( [16] ) is taken over all components i for which the resonance condition (16) is 
satisfied (for a given order p). The "source term" 3> pq on the right-hand side of Eq. (46) is 
given by the expression 

r 

* P q = f P q - h p'q' X^ u rg(p-p'+<5 3 )(q-q')- ( 47 ) 
p'q' i =1 

in terms of f pq and lower order coefficients h P ' q ' and g P " q " , given by 

u*-g P " q " = | Uj *' fp " q " if P"' A = A ^ (48) 
3 1 otherwise. 

The vectors f pq are coefficients in an expansion of the vector field f (x, /x) on the unfolded 
center manifold W c (^), 

f(y + h(y,/i),M) = X;W P /* q - (49) 
pq 
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In practice, this expansion is obtained by substituting the expansion ( |44| ) in a Taylor expan- 
sion of f (x, ft) in terms of the usual multilinear forms, e.g. 

F X x(x,x)- 2^ q g XiXj, F xxx (x,x,x)_ ^ d Xi d Xj dx k XiXjXk ' 

1 i,j.k—l 

50 

s <9F " s <9 2 F 

i=l ^ i=l j=l ^ 



Any coefficient f pq of the expansion (|4£|) is then obtained as an expression that may 
involve coefficients h p / q /. The important point now is that only lower orders (p',q') appear 
for lw 

in fpq- So it is always possible to solve Eqs. (|45|) — (ff7|) order by order. 
As a result, the amplitude transformation ( ffi| ) is obtained explicitly. This procedure 
provides a particular representation of the unfolded center manifold W c (fj.) in terms of the 
amplitudes j/j. Some simple examples of the calculation of the coefficients f pq and explicit 
solutions of the linear equation ( f45| ) can be found in |Q| . 

B Proof of Eqs. (03) and ([13]) 

We use the expansion in Eq. (|4i| ) for h(y,/x) together with Eq. (|l4|). For the left-hand side 
of Eq. , we then get 

£^h pqZ P^]T^i fe^-^dO. (51) 

pq j=l Z 3 T J 

Here a term of order (p,q) vanishes because of the integral over 9 unless p-A = A^. So the 
sum contains resonant terms only, but each of these vanishes because resonant components 
of h pq vanish, u*h pq = 0, as Eq. ( 53 ) and Eq. ( fl5| ) show. This demonstrates the identity 
in Eq. @. 



For Eq. (|13[), we use Eq. ( |44] ) and first observe that 

J-y-|=J-o Je -z-|(c-.z) = (52) 

From the expansion ( fill ) of the transformation h(y, /z), we therefore get from the left-hand 
side of Eq. (|l|) 

£ <-(J - p - A) •hpqzP/Zii f e^ x -^ d9. (53) 

Here only resonant terms contribute to the sum because of the integral, but these also vanish 
since u*-(J — p-A) = 0. The whole sum is therefore identically zero, so Eq. (13) is proved. 



C Relation to the CGLE 

Comparsion of the DSHE (E7f with the CGLE Eq. (M) contains a subtlety, which we now 



discuss. Equation (27a) has the same form as Eq. (]26[) apart from the coupling term g\wz. 
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But the coefficients g\ and g also differ. They are defined in Tables [j] and ^ by very similar 
expressions, namely 

gi = u*-F xx (u,hnoo) + u*-F xx (u,h 20 oo) + ^u* F xxx (u, u, u), (54) 
g = u*-F xx (u,hno) + u*-F xx (u,h 20 o) + |u*-F xxx (u,u,u), (55) 

which differ only by hno replacing hnoo and h2oo replacing li2ooo- The latter pair are 
identical because they are unique solutions to identical linear equations. But hno 7^ hnoo 
because (p, q) = (1, 1, 0, 0) is resonant for the real mode in the fold-Hopf bifurcation, whereas 
(p, q) = (1,1,0) is nonresonant in the Hopf bifurcation. So the two vector coefficients are 
defined through two slightly different linear equations: 

-J-hnoo = F xx (u,u) - c v, (56) 
-J-hno =F xx (u,u). (57) 



Here, Co is the coefficient of \w\ 2 in Eq. (27b). The second term of Eq. ( |5q ) arises through 
the projection operator Q3 in Table [j] 

The physical significance is the following. In the CGLE, effects of components of the vector 
field F out of the plane of oscillation can only be taken into account through the amplitude 
transformation. In the DSHE, the dynamics of the slow mode is directly accounted for in 
the amplitude equation. Consequently, it does not appear in the amplitude transformation, 
and as a result, the coefficients g\ and g are in fact quite different. Compare their numerical 
values for the chosen operating point of the 4D Oregonator exhibited in Tables || and |(| Here, 
it is interesting to note that using g\ in place of g in the CGLE would place the equation 
well inside the region of absolute instability (see Fig. |l|d). 

However, this comparison of parameters is misleading, since it neglects the effect of the 
term gowz. To illuminate the issue, we first find an explicit relation between hnoo and hno 
from the difference between Eqs. ( p6| ) and (p7|), 

J • (hnoo ~ hno) = c v. (58) 

Since J is nonsingular , the solution is 

hnoo - hno = c J _1 -v = -^-v. (59) 

Ao 

We can then use this result to derive an explicit relation between g\ and g. Solving Eq. (|59| ) 
for hnoo and substituting the result in Eq. (jHJ), we get, using the bilinearity of F xx (u, hnoo) 
and the identity I12000 = h2oo: 

.gi =.g + ^u*-F xx (u,v) =g+^g . (60) 

Ao Ao 

We can get an indication of the effect of the term g w z by examining the artificial situation 
where z is constrained to be homogeneous and stationary, so 

z = -fH 2 (61) 

Ao 



from Eq. (27b). Substitution of this result in Eq. (E7aj) then yields 



w = (Tifiw + (gi - — )\w\ 2 w + d w V 2 w, (62) 
Ao 
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which reduces to the CGLE, exactly, when Eq. ( |60| ) and the identity d w = d are employed. 
Thus it is the variation of z with space and time — the dynamics of the slow mode — that 
causes the deviation from CGLE behavior, which of course is no surprise. 

In practice, we regain CGLE behavior if Ao becomes sufficiently large and negative im- 
plying that g\ and g almost coalesce (according to Eq. (|6fj|)) and the linear term Xqz in fl27b| ) 
ensures that \z\ and hence |go W2; | never becomes large. The magnitude of Ao should be 
compared with ^Re{<7i} first of all, and a condition for using the CGLE is |Aq| 3> ^Re{cri}. 
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Captions of tables and figures 

Figure |j: (a) : Bifurcation diagram showing the locations of supercritical Hopf bifurcations in 
the 4D Oregonator model in the plane of the two parameters [BrO^] and [H + ]. (b): Variation 
of the dimensionless diffusion parameter a along the bifurcation curve in (a), plotted as a 
function of the parameter [BrO^~]. (c): Variation of the dimensionless parameter /3 along the 
branches in (a) plotted as a function of the parameter [BrO^]. (d): Projections of the curves 
in (b) and (c) onto the (a,/3)-plane together with the curves corresponding to the Eckhaus 
(EH) and absolute instability (Al) curves respectively. The Benjamin-Feir instability line 
(BF) is also shown. 

Figure |2|: Solution to the 4D Oregonator PDE at a finite distance from the bifurcation 
corresponding to an amplitude, which is 7.5% of the value of the stationary Ce 4+ concentra- 
tion. The three rows show the amplitude \w\ and two selected concentrations at the times 
shown below each column. Integrations were made on a 256x256 grid with no-flux boundary 
conditions and physical dimensions 4 cm in both directions. 

Figure |: Solution to the CGLE corresponding to the one for the 4D Oregonator shown 
in Fig. |, The number below each column shows the time elapsed for that particular state. 
Integrations were made on a 256x256 grid with no-flux boundary conditions and physical 
dimensions 4 cm in both directions. The relative variation of [BrMA] is extremely small 
(compare with Table 0) . It is represented correctly by the colors but cannot be distinguished 
on the corresponding colorscale. 

Figure Solution to the 4D Oregonator PDE at a finite distance from the bifurcation point 
corresponding to an amplitude, 7.5% of the value of the stationary Ce 4+ concentration. The 
number below each column shows the time elapsed for that particular state. The solution 
is the same as in Fig. ^, but the amplitude z is exhibited instead of \w\. Integrations were 
made on a 256x256 grid with no-flux boundary conditions and physical dimensions 4 cm in 
both directions. 

Figure |B[ Solution to the DSHE corresponding to that for the 4D Oregonator shown in 
Fig. [|. The number below each column shows the time elapsed for that particular state. 
Integrations were made on a 256x256 grid with no-flux boundary conditions and physical 
dimensions 4 cm in both directions. 

Figure ^: Time series showing the local variation of Re w (top) and z (bottom) at a particular 
grid point in the course of the pattern development shown in Fig. |[ Observe that the time 
variation of z is much slower than for flew. 

Figure fr|: (a): Solution to Eq. (42) approaching the stable focus (r s ,z s ). 

Table [jj: Formulas for calculating the coefficients of the amplitude transformation and am- 
plitude equation for the fold-Hopf bifurcation. At the bifurcation, the Jacobian J has three 
eigenvectors u, u, and v and left eigenvectors u* , u* , and v* corresponding to the three crit- 
ical eigenvalues Ai = A2 = iwo and A3 = 0. The amplitude transformation x = y + h(y, /x), 
y = uwe 1Uot + c.c. + vz, transforms a solution w(r,t), z(r,t) of the amplitude equation to 
the motion x(r, t) on the unfolded center manifold for the dynamical system. The vector 
coefficients h pqrs are determined as solutions to the linear equations indicated, in terms of 
the derivatives of the vector field F. The coefficients of the amplitude equation can then 
be found through the explicit expressions indicated, in terms of the derivatives of F and 
h pqrs . For any x e K n , the projections Qi and Q3 are defined as Qi-x = x (u*-x)u and 
Q,3- x = x (v*-x)v respectively. 
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Table ||: Formulas for calculating the coefficients of the amplitude transformation and am- 
plitude equation for the Hopf bifurcation. At the bifurcation, the Jacobian J has two complex 
conjugate eigenvectors u and u and left eigenvectors u* and u* corresponding to critical eigen- 
values Ai = A2 = iu)Q. The amplitude transformation x = y + h(y, /x), y = uwe lu>at + c.c, 
transforms a solution w(r, t) of the amplitude equation to the motion x(r, t) on the unfolded 
center manifold for the dynamical system. The vector coefficients h pqs are determined as 
solutions to the linear equations indicated, in terms of the derivatives of the vector field F. 
The coefficients of the amplitude equation can then be found through the explicit expressions 
indicated, in terms of the derivatives of F and h pqs . For any xgl", the projcetion Qi is 
defined as Qi x = x (u* -x)u. 

Table |^: Rate constants fei, . . . ,kr, constant malonic acid concentration M, and diffusion 
constants D x , D Y , D z , and D v used for numerical simulations for the 4D Oregonator 
reaction-diffusion equation (j3|) . The values of the rate constants are based on || but have 
been modified slightly in the calculations presented here. Diffusion constants Dx, ■ ■ ■ ,Dz 
are from , whereas the value of Djj was estimated by an interpolation scheme described 
in H p. 196]. 

Table |^: Hopf bifurcation point for the 4D Oregonator used as a reference for the operating 
point actually used in the calculations. The table shows the parameters, the stationary point, 
and the eigenvalues ±icJo, A3, and A4 of the Jacobian matrix at that point. 

Table ^ Ginzburg-Landau parameters at the Hopf bifurcation, Table for the 4D Orego- 
nator Eq. (|35|) with parameters shown in Table @. 



Table Ej|: Coefficients for the distributed slow-Hopf equation ( |27D at the Hopf bifurcation, 
Table [5 , for the 4D Oregonator Eq. ( |35| ) with parameters shown in Table ||. 

Table @: Left and right eigenvectors of the Jacobian matrix at the Hopf bifurcation, Table ||, 
for the 4D Oregonator Eq. (35) with parameters shown in Table ||. The eigenvectors u* (left) 
and v (right) correspond to eigenvalues iwo and A3 respectively. They illuminate the special 
role played by BrMA in the 4D Oregonator. 




Figure 1: 
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Figure 6: 




Figure 7: 
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Fold-Ifopf Bifurcation 


Transformation 

x = y + h(y.^) 


x = uip + uif> + wz + h 2 ooo</" 2 + hnoo|-0| 2 + h 200>/' 2 + h o2oz 2 + hioioV'z + h no'/' z + 
h 30 oo>/> 3 + h 2 ioo l"0| 2 ^ + h 12 oo\i>\ 2 i> + ho30oV" 3 + h o 3 oz 3 + h 2 owif> 2 z + huw\i>\ 2 z + 
h 02 wip 2 z + h 102 oi>z 2 + h 012 oi>z 2 + (hnnni + h WO iip + hnioiV' + hooliz)^, i> = we'" ot 




Second order 


Linear equations for h pqrs 


w 2 


-(J - 2 iwrjl) • h 2 ooo = |F XX (u. u) h 02 oo = h 2 ooo 


M 2 


-J-hnoo = -Q 3 -F xx (u,u) 


wz 


-(J - iu u I)-hioio = Qi-F xx (u, v) u*-hoo20 = hono = hioio 


z 2 


-J-hnn20 = 2"Q3-Fx x (v, v) V*-h 020 =0 




Third order 


Linear equations for h pqrs 


w 3 


-(J - 3iw I)-h 3 ooo = F xx (u, h 2 ooo) + gF xxx (u, u, u) ho300 = h 30 oo 


\w\ 2 w 


-(J - iwoI)-h2ioo = Qi-(F xx (u,hnoo) + _ 

U* -112100 = hi200 = h 2 ioo 

F xx (u, h 2 ooo) + |F xxx (u, u.u)) 


wz 2 


- (J - i w I) • hio20 = Qi ■ (F xx (u, h 020 ) + _ 

hoi20 = hi020 

F xx (v, hioio) + |F xxx (u, v, v)) 


\w\ 2 z 


-J-hmo = Q 3 -(F xx (u, h iio) + F xx (u, hioio) + F xx (v, hno ) + F xxx (u, u, v)) 


W 2 Z 


-(J - 2iu I)-h 2 oio = F xx (u, hioio) + F xx (v, h 2 ooo) + _ 

ho2io = h 2 oio 

j-Fxxxfu, u, v) - g 2 h2ooo 


z 3 


-J-h 030 = Q3-(F xx (v,h 00 2o) + |F xxx (v, v, v)) v*'h o 3 o =0 




Unfoldings 


Linear equations for lipqrs 


M 


— J-hoooi = Q 3 -F^ v*-hoooi =0 


fllO 


-(J - iwoI)-hiooi = Qi'(F XfI -u + F xx (u, hoooi)) u*-hiooi=0 h O i i = hiooi 


flZ 


-J-hoon = Q 3 -(F X)J -v + F xx (v,h ooi)) v*-h on =0 




Amplitude 
equation 


w = ainw + gowz + gi \w\ 2 w + g^wz 2 + d w ^7 2 w, 

d w = u* ■ D ■ u, d z = v* ■ D ■ v 

z = pofi + cq\w 2 + c±z 2 + C2 \w\ 2 z + c 3 js 3 -f d z V 2 z. 


Resonant 
coefficients 


go = u* F xx (u, v) 

gi = u* F xx (u, hnoo) + u* F xx (u, h 2 ooo) + Ju* F xxx (u, u, u) 
g 2 = u* F xx (u, h o2o) + "* -F xx (v, hioio) + |u*-F xxx (u, v, v) 
c = v*-F xx (u,u), ci = |v*-F xx (v,v) 

c 2 = v*F xx (u, hono) + v* -F xx (u, hioio) + v*-F xx (v, hnoo) + v* F xxx (u, u, v) 

c 3 = v*-F xx (v,h 00 2o) + gv* -F xxx (v, v, v) 

cri = u* F XfI u + u* F xx (u, hoooi), Pa =v*-F tl 



Table 1: 
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Hopf Bifurcation 


Transformation 

x = y + hfy, fjt) 


x = uV> + utf> + h2oo'0 2 + h 1 1 o 1 1 2 + ho20ip 2 + h3ooV" 3 + h 210 \il> 2 i/> + 

hi2o|'i/'| 2 i/J + ho30'/' 3 + hnoiM + (hioi - !/ 1 + hoiii/')/*. i> = we 1 "' ' 




Second order 


Linear equations for h pqs 


w 


-(J - 2iwoI)-h200 = jF X x(u,u) h 20 = h 2 oo 


\w\ 2 


-J-hno = F xx (u,u) 




third order 


Linear equations for h P q S 


w 3 


-(J - 3iwoI)-h 30 o = F xx (u, h 2 oo) + §F xxx (u, u, u) h 030 = h 30 o 


\w\ 2 w 


-(J - iw I)'h2io = Qi -(F xx (u, huo) + _ 

u*-h 2 io = hi20 = h 2 io 

F xx (u,h200) + 5 F x xx (u,U, u)) 




Unfoldings 


Linear equations for hpq S 


M 


-J-hnoi = F,, 




-(J - iw I)-hioi = Qi -(F xfI u + F xx (u, hooi)) u*-hioi=0 h ii = h 10 i 




Amplitude 
equation 


w — aifiw + g\w\^w -f a"V 2 w d = u*-D-u 


Resonant 
coefficients 


g = u*F xx (u, huo) + u* F xx (u, I1200) + fu* -F xxx (u, u, u) 
a i = u*F xfl u + u*F xx (u, hooi) 



Table 2: 



Constant 


Value 




Constant 


Value 


fci/M'^s -1 


1.6 




fc 7 /M _i 


0.0003 




fc 2 /M _2 s _1 


2.5 x 10 6 




M/M 


0.44 




fc 3 /M _2 s _1 


33.0 




Djc/cm 2 s _1 


1.0 x 10" 


-5 


fc 4 /M- 1 s- 1 


3.0 x 10 3 




Dy 1 'cm 2 s _1 


1.6 x 10" 


-5 


fc 5 /M- 1 s- 1 


30.0 






0.6 x 10~ 


-5 


fce/M^s" 1 


0.18 




D v /cm 2 *- 1 


6.7 x 10" 


-6 



Table 3: 



Reaction-Diffusion Systems with Slow Real Modes 



Parameter 

[BrO^Hopf/M 

[H + ] H o P f/M 

[HBrO] /M 

[Br-] /M 

[Ce 4 +] /M 

[BrMA] /M 

A3/S- 1 
A4/S- 1 



Value 
03662 
1.3416 

9.8679 xlO" 7 
7.0929 xlCT 6 
1.3185 xl0~ 5 
7.8262 xl0~ 2 
2.1911 
-3.0731 xlO~ 4 
-10.4379 



Table 4: 



Parameter 




Value 




-2.57 


+ i 1.14 


g/lO^M^s^ 1 


2.19 


+ i2.09 


d/lO^cm^" 1 


1.00 


- i 2.10 


a 


0.96 




P 


-0.21 





Table 5: 



Parameter 


Value 


Ao/s" 1 


-3.07 x 10~ 4 




-2.57 + i 1.14 


g /10 5 M- 1 s- 1 


0.64 + il.00 


g 1 /lO w M~' 2 s~ 1 


1.09 + i3.80 


co/10 M-^- 1 


-5.25 


d w /lO- 5 cm 2 s- 1 


1.00 - i 2.10 


da/lO -6 cm 2 s _1 


6.70 



Table 6: 



Species 


Reu* 


Imu* 


V 


HB1O2 


7.546 x lO" 1 


-1.843 


3.382 x lO" 3 


Br" 


1.699 x 10" 1 


8.390 x 10" 1 


-1.138 x 10~ 2 


Ce 4 + 


4.665 x lO" 1 


2.519 x 10" 1 


1.0 


BrMA 


7.975 x 10" 5 


4.032 x 10~ 5 


-5.857 x 10 3 



Table 7: 



